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A simple coupled oscillator system with chemotaxis is introduced to study morphogenesis of 
cellular slime molds. The model successfuly explains the migration of pseudoplasmodium which 
has been experimentally predicted to be lead by cells with higher intrinsic frequencies. Results 
obtained predict that its velocity attains its maximum value in the interface region between total 
locking and partial locking and also suggest possible roles played by partial synchrony during 
multicellular development. 
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§1. Introduction 

Inspired by the collective behavior and rhythmicity of 
biological systems, synchronization of coupled limit-cycle 
oscillators with a frequency distribution has been studied 
using.a.s.KStem with all-to-all interaction in the following 
forml» 
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Wj), (1.1) 



where tUj is an intrinsic frequency of the oscillator j, Wj 
is a complex variable and (') = 4rQ- The system has 
served well for theoretical studies on frequency entrain- 
ment and critical behaviors in a system far from equi- 
librium. Despite its original aim, however, very little 
has been discussed on the system's applicability and a 
connection to biological systems. 

Development of the cellular slime mold Dictyostelium 
discoideum provides us with an ideal example that needs 
to be investigated in this respect, since it's chemoat- 
tractant secretion exhibits limit-cy.cle oscillations in the 
vicinity of a Hopf bifurcation points' wittua difference in 
intrinsic frequencies and synchronization© 

These eukaryotic cells feed on bacteria and grow by bi- 
nary fission. Deprived of food, they initiate aggregation 
by emitting cAMP as a chemoattractant while simultane- 
ously making a directed motion against its surrounding 
gradient. Each aggregation territory consists of 10 3 to 
10 5 amoebae that together form a spherical mound whose 
nipple like apex zone consists of differentiated prestalk 
cells. The tipped mound elongates vertically to form a 
slug like pseudoplasmodium which migrates to.a suitable 
environment where it forms a fruiting body.BLP 

To study a self-organizing process of cells in such sys- 
tems, one must understand the nature of collective be- 
havior in a population of motile oscillators. In the follow- 
ing, we first int rod uce a system that incorporates chemo- 
taxis into eq. (IT) and present the system's overall be- 



havior. Then we discuss our results and their implica- 
tions on development of a cellular slime mold. 

§2. Equations 

The system is derived from a linear diffusion equation 
for two chemical species denoted by a complex variable 
Z(r, t) and equations for chemotaxis of cell j at r = rj. 
They are 



d t Z(r,t) = Z?V 2 Z(r,i), 
j-j(t) = aVReZ(r,i), 



(2.1) 
(2.2) 



where D is a diffusion constant, a is a chemotactic sen- 
sitivity coefficient. 

Cell j (= 1,2,- ■ -,N) is represented by a region 
[0 < |r — rj| < ro] on which a boundary condition 
representing the metabolism T of intracellular chemi- 
cal species in a complex variable Wj (t) is imposed. The 
boundary conditions are expressed as 

lim Z(r,t) = 0, (2.3) 

r— »oo 

Z(rj+r ,t) = Wj(t)+C, (2.4) 

Wj(t)=F(Wj)+iP({W k }). (2.5) 

Here ip({Wk}) is an abbreviation for the interaction term 
ip({Wi,W2, • ■ -,Wn}), and C is a constant background 
level of Z(r,t). Solving cq. (2.1) under N boundaries 
(2.4) would yield a field Z(r,t) in an integral equation 



Z(r,t) 
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*(|r(t) - r fc (r)|,t - r)(W fe (r) + C)dr, 

(2.6) 

where the diffusion kernel $ is either a Gaussian or 
Bessel-type function according to the dimension of r. 

In order to simplify the system, let us assume $(|r(t) — 
rk(r)\,t — t) — > 5(t — t) in the limit of \r — r^j — > 0. 
Assuming we had T that yields the Hopf bifurcation 
normal form under a mean-field coupling ip({Wk}) oc 
(Z(rj,t) - Z(rj + r ,t)); eqs.(2.2) and (2.5) could be 
rewritten as 
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where uij and A independently determine the intrinsic frequency and the amplitude of the oscillation. In a weakly 
coupled regime, we could neglect the amplitude effect, therefore allowing eqs.(2.7) and (2.8) to be further reduced 
to a coupled phase model in the following form, 

e N ft 

M*) = ^ + J^J E / $ d r i(*) -r k (r)\,t- T ) sm(<Mr) - 0j(*))dr, (2.9) 
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where c is the real part of C. In the transformation from 
eq. (2.7) to eq. (2.9), we have fixed A = 1. In addition, 
the deceleration term and a periodic change of chcmjQ- 
taxis sensitivity were added to our previous model© 
Here, the sensitivity function 7 and the short range in- 
teraction Vd are defined by 



7(0) = a(l 



(2.11) 



Vd = - 
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(2.12) 
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where (3 and k are both positive constants. 

By imposing a boundary condition on each cell that 
is in the form of an ordinary differential equation, the 
spatially extended system in eqs. (2.1) and (2.2) is now 
reduced to a set of integro-differential equations that 
tracks the time development of those cell boundaries. 
Notice that eq. (2.9) provides us with a general scheme 
that incorporates, spatial dependency in the well-studied 
phase model.li^'tLy) This is a novel approach to reaction- 
diffusion systems applicable to those that exhibit nonlin- 
earity localized on boundaries. 

§3. Method 

For the sake of numerical analysis, the kernel is sim- 
plified down to a point where it is not distinguishable 
from an one dimensional kernel except by the factor of 
y/^f and a shift in the origin of the exponent. Further- 
more, it is multiplied by a step function which roughly 
incorporates the effect of degradation of the chemoat- 
tractant by the enzyme in the extracellular substratum. 
Therefore At could be considered as an average life span 
of the chemoattractant. The precise form used in the 
calculations is as follows: 



$(r,i- t) 
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where 9(i) denotes a Heaviside step function. 

Numerical studies on eqs. (2.9) and (2.10) were per- 
formed using the fourth-order Rnnge-Kutta method in 
addition to the semiopen formulaEf for the integral terms. 
Parameters D, ro and n were fixed at unity through- 
out latter calculations. Other parameters are c = 2.0, 
LJj = 1.0 + (j - l)Aw, (3 = 0.01, At = 2.0 unless stated 
otherwise. 

We have employed constant initial values for the inter- 
val of [—At, 0], and all results were obtained from a fixed 
step size of h=0.01. Some calculations where N = 2 were 
checked for accuracy using h=0.001 and h=0.0001. 

§4. Migration in TV = 2 

We will first describe the simplest case where N = 
2 to characterize attractors in the system. In order to 
carry out a steady-state approximation, chemotaxis will 
be confined to one dimension, and the adaptation will 
also not be considered (k = 0). 

Suppose oscillators were entrained with a constant 
phase difference tjj — <p\ — <f>2- When X2 — x\ = 2tq = 
const., the equation describing the position of a centroid 
x c = (xi + X2)j2 could be approximated by 



[ V$(2r ,T / )sin 9 !)i(<-T / )dr', 
Jo 



(4.1) 



where we have neglected the second- and higher-order 
terms of i[> and had a change in variable from r to t' = 
t — t. The mean cluster velocity v c =< x c > therefore 
would be proportional to ip. 

We see that v c is also proportional to Aoj from the 
locking condition -0 = 0. Applying the same approxima- 
tion as in eq. (4.1), 

i> ~ -Auj - 2i>e [ $(2ro,T')cosw*r'dT', (4.2) 



where uj* 
locked, 



is the entrained frequency. When the phase is 
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Fig. 1. and the cluster velocity v c (ui = 1.0, u>2 = uii + 

Aw,e = 1.0 and a = 1.0). v c was obtained using the method of 
least squares. 



Fig. 3. e and the order parameter R (a = 1.0 and Alu = 0.01). 
The minimum is plotted to withdraw fluctuations due to small 
N. 
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Fig. 2. a and the cluster velocity for different e (uii = 1.0,^2 = 
1.1 and a = 1.0). The initial condition was chosen so that 
xi (0) < x 2 {0). 



must be satisfied. From eqs.(4.1) and (4.3), one obtains 
|u c | oc aAuj/e. 

The parameter dependencies given here were con- 
firmed by numerical simulations for N — 2 with k = 1.0. 
For a sufficiently large a, the oscillator with larger u>j is 
advanced in phase and leads migration as one can see in 
Figs. 1 and 2. Not only a but Aui also increases the clus- 
ter velocity v c . Note that there is no net migration of a 
cluster when Aui = 0. Figure 2 also indicates that there 
is velocity proportional to a (but not to 1 /e) even when 
frequency is not entrained (e = 0.1) due to the deceler- 
ation term which was not present- in Aizawa-Kohyama 
(A-K) model discussed previously M> 



§5. Numerical Results for N = 20 

In order to characterize coherence in the phase dy- 
namics, an order parameter R — jj J2f=i i s plotted 
against e in Fig. 3. Since we are dealing with a very small 
number of N, \R\ does not approach zero as e — ► 0. Note 



Fig. 4. e and cluster velocity (a = 1.0 and Au; = 0.01). 



that in Fig. 3, the minimum of \R\ was plotted instead 
of, for example, < \R\ > (<> denotes time averaging). 

A cluster shows a directed migration when oscillators 
are entrained into a common frequency. The instanta- 



neous velocity of a centroid defined by r c 
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is plotted against e in Fig. 4. The apparent discontinuity 
at the onset of total entrainment suggests that the clus- 
ter velocity, too, may be taken as an order parameter. 



The relation between the polarity of a cluster and the 
migration direction could be understood by introducing 
Z(6) which is defined by 



N 



Z{6) =Y,km k {6), 



(5.1) 
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where m^s are natural numbers {1,2,3,- • -,N} repre- 
senting the distance of oscillator j from a given point 
outside the cluster (x,y) = (r cos 9, r sin 9); 1 is assigned 
to the closest and 20 to the furthest oscillator. From the 
above definition, Z{&) takes the maximum value in the 
direction 9 (measured from the centroid) where oscilla- 
tors with a smaller u)j are positioned. 
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one exhibited, by subpopulation of cells at the onset of 
culmination^ could be realized without any secondary 
chemoattractant. It would be interesting to see whether 
cells make use of such effect caused by different coupling 
ranges. This could be controlled either by the extra- 
cellular enzyme concentration or by the surface receptor 
occupancy. 

The relation between the coupling strength e and the 
instantaneous cluster velocity Jhas partly confirmed the 
result obtained by A-K mode© with its peak centered 
around the critical region. Below e = 7.5, v c continues 
to fluctuate with its direction not fixed which prevents 
us from obtaining its mean. 

We have observed that in a partially synchronized 
state, locked oscillators emerge specifically in the ante- 
rior section of a cluster (data not shown) which implies 
that desynchronization plays a possibile role in cell-type 
differentiation. In general, / _At fluctuations are to be 
observed in such interface region. A detailed analysis on 
such fluctuation and its effect will be provided elsewhere. 

Due to the coupling and chemotaxis in the form of 
integral equation, we have confined our present report on 
results obtained in the case of N = 20 with fixed tq and 
D. The diffusion characteristic makes it difficult to vary 
the ratio r^/D since decreasing it requires smaller step 
size h and increasing it requires larger At. Future work 
will also address the construction of a simpler scheme 
that makes the analysis more feasible. 



Fig. 5. Cluster polarity and migration direction (e = 30.0, a = 
1.0, Au; = 0.01). Vertical lines indicate the direction of cluster 
migration. Other parameters are At = 2.0 (a) and At = 6.0 (b). 
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Figure 5 displays Z(8) and the direction of cluster mi- 
gration. When e and At are small but large enough to 
synchronize oscillators, a cluster heads toward the direc- 
tion where more oscillators with smaller uij are located. 
The orientation reverses as both e and At are increased 
so as to make the coupling long range. 

§6. Discussion 

A simple coupled oscillator model of cell aggregates 
was derived from a linear diffusion equation with time- 
dependent boundaries. The approximation in N = 2 
and numerical analysis carried out for N — 20 revealed 
that synchrony in such a population of oscillators with 
chemotaxis results in migration as a whole. In addition 
to the migration direction that agrees with an experi- 
mental prediction, there are some implications from our 
present work concerning the possible order parameter. 

We showed in Fig. 5 that in the case of N = 20, oscil- 
lators with larger ujj also lead the cluster translocation 
if At is sufficiently large. In-the light of the prediction 
from suspension experiments© that cells with higher fre- 
quency constitute the anterior of a slug, it may be con- 
cluded that cell to cell interaction is not local but rather 
long range. The opposite migratory direction predicted 
by a small At suggests that a reverse flow such as the 
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